####################################################################################
#
# Replication file for: Ingram, Matthew C. 2016. "Mandates, Geography, and Networks: 
#     Diffusion of Criminal Procedure Reform in Mexico"
#     Latin American Politics and Society 58: 121-145.
#     DOI: http://dx.doi.org/10.1111/j.1548-2456.2016.00301.x
# Author: Matthew C. Ingram
# Affiliation: University at Albany, SUNY
# contact: mingram@albany.edu
# NOTE: this file complements other command files (see replication materials)
# 
####################################################################################

####################################################################################
# SET ENVIRONMENT: packages and directory structure
####################################################################################

# if needed, install the following packages prior to opening libraries
# e.g., install.packages('spdep')
library(spdep)
library(maptools)
library(rgdal)
library(spgwr)
library(RColorBrewer)
library(foreign)

# set working directory (note: this is my filepath; change to your own working directory)
path <- "C:/Users/MI122167/Dropbox/Criminal Procedure/pubs/Ingram2016LAPS/replicationmats"
setwd(path)

# check to see if sub-directories exist; create them if they are not there already
dir.create(paste(path,"/","data", sep=""))
dir.create(paste(path,"/","figures", sep=""))

####################################################################################
# read data
####################################################################################

# load Mexico state shapefile
mex <- readOGR(dsn="shapefiles",layer="Mexico_congruence")

mexdata <- mex@data

# Load spatial weights data; could also generate these in R, but already had these on file

mexnbq1 <- read.gal("weights/LASA2014_q1.gal", region.id=mexdata$ST_ID)
# check object to confirm
mexnbq1
summary(mexnbq1)
head(mexnbq1)

# ensure it is symmetric
mexnbq1sym <- make.sym.nb(mexnbq1)
mexnbq1sym

mexnbnet <- read.gal("weights/LASA2014_partynet2010.gal", region.id=mexdata$ST_ID)
mexnbnetsym <- make.sym.nb(mexnbnet)
mexnbnetsym

# New DVs and lags
w <- nb2mat(mexnbq1sym)
dim(w) # shoudl return 32 x 32
w
str(w)


# using library(foreign)
data <- read.dta("data/2newdvs2015.dta")
head(data)

wi4l02 <- w%*%data$i4l_02
wi4l03 <- w%*%data$i4l_03
wi4l04 <- w%*%data$i4l_04
wi4l05 <- w%*%data$i4l_05
wi4l06 <- w%*%data$i4l_06
wi4l07 <- w%*%data$i4l_07
wi4l08 <- w%*%data$i4l_08
wi4l09 <- w%*%data$i4l_09
wi4l10 <- w%*%data$i4l_10
wi4l11 <- w%*%data$i4l_11

wi4i02 <- w%*%data$i4i_02
wi4i03 <- w%*%data$i4i_03
wi4i04 <- w%*%data$i4i_04
wi4i05 <- w%*%data$i4i_05
wi4i06 <- w%*%data$i4i_06
wi4i07 <- w%*%data$i4i_07
wi4i08 <- w%*%data$i4i_08
wi4i09 <- w%*%data$i4i_09
wi4i10 <- w%*%data$i4i_10
wi4i11 <- w%*%data$i4i_11

wi4lwide <- (cbind(wi4l02,wi4l03,wi4l04,wi4l05,wi4l06,wi4l07,wi4l08,wi4l09,wi4l10,wi4l11))
wi4lwide
dim(wi4lwide)
str(wi4lwide)
# convert dgeMatrix to regular matrix and then to data frame
wi4lwide <- as.data.frame(as.matrix(wi4lwide))

wi4ltall <- reshape(wi4lwide,
                    varying = c("V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8", "V9", "V10"), 
                    v.names = "q1i4l",
                    timevar = "year", 
                    times = c("2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011"), 
                    new.row.names = 1:(32*10),
                    direction = "long")
wi4ltall
summary(wi4ltall[,"q1i4l"])

wi4iwide <- (cbind(wi4i02,wi4i03,wi4i04,wi4i05,wi4i06,wi4i07,wi4i08,wi4i09,wi4i10,wi4i11))
wi4iwide
dim(wi4iwide)
str(wi4iwide)
# convert dgeMatrix to regular matrix and then to data frame
wi4iwide <- as.data.frame(as.matrix(wi4iwide))

wi4itall <- reshape(wi4iwide,
                     varying = c("V1", "V2", "V3", "V4", "V5", "V6", "V7", "V8", "V9", "V10"), 
                     v.names = "q1i4i",
                     timevar = "year", 
                     times = c("2002", "2003", "2004", "2005", "2006", "2007", "2008", "2009", "2010", "2011"), 
                     new.row.names = 1:(32*10),
                     direction = "long")
wi4itall
summary(wi4itall[,"q1i4i"])

# col bind two Queen-1 lags of i4l and i4i
q1lags2 <- cbind(wi4ltall, wi4itall["q1i4i"])
names(q1lags2)[names(q1lags2)=="id"] <- "code"
head(q1lags2)

# write out new data to file
write.dta(q1lags2, "data/2newds2015_q1lags.dta")


#######################################################################
# Figure 1 (May 2015)
#######################################################################

mex2 <- readOGR(dsn="shapefiles", "LASA2014_mergedDVs")

# following code generates Figure 1 and saves it to file as 300dpi .tif file
tiff(file="figures/figure1.tif", width=6, height=6, units="in", res=300)
spplot(mex2, "i42011", main="2011 Reform Index", col.regions=brewer.pal(7, "Blues"), cuts=6,scales=list(draw = FALSE))
dev.off()


# load party-specific shapefiles and coordinates
mexpan <- readOGR(dsn="shapefiles", "Mexico_congruence_PAN2010")
mexprd <- readOGR(dsn="shapefiles", "Mexico_congruence_PRD2010")
mexpri <- readOGR(dsn="shapefiles", "Mexico_congruence_PRI2010")

# load party-specific weights
mexnbnet2010pan <- read.gal("weights/LASA2014_partynet2010_panonly.gal", override.id=T)
mexnbnet2010prd <- read.gal("weights/LASA2014_partynet2010_prdonly.gal", override.id=T)
mexnbnet2010pri <- read.gal("weights/LASA2014_partynet2010_prionly.gal", override.id=T)

#######################################################################
# Figures 3 and 4 (May 2015)
#######################################################################

# simple version of Figure 3 plot with queen ties
coords <- coordinates(mex2)

png("figures/figure3.png")
par(mar=c(0.1,0,1,0), omi=c(0,0,0,0))
plot(mex2, border="grey", axes=F)
plot(mexnbq1, coords, col="red", pch=".", add=T)
text(coords, label=mex$ST_ID, cex=.75)
title("Spatial Proximity (Queen-1)")
dev.off()

# Figure 3; code to generate final 300dpi .tif file

tiff(file="figures/figure3.tif", width=6, height=6, units="in", res=300)
par(mar=c(0.1,0,1,0), omi=c(0,0,0,0))
plot(mex2, border="grey", axes=F)
plot(mexnbq1, coords, col="red", pch=".", add=T)
text(coords, label=mex$ST_ID, cex=.75)
title("Spatial Proximity (Queen-1)")
par(mar=c(5, 4, 4, 2)+0.1)
dev.off()


# Figure 4; preparation

coordspan <- coordinates(mexpan)
coordsprd <- coordinates(mexprd)
coordspri <- coordinates(mexpri)

# Figure 4; code to generate final figure and save 300dpi .tif file

tiff("figures/figure4.tif", width=6, height=6, units="in", res=300)
par(mar=c(0,0,1,0), omi=c(0,0,0,0)) # margins (bottom, left, top, right) for each graphs
layout(matrix(c(1,2,3,4,4,4), 2, 3, byrow = TRUE), #"2" establishes number of rows, "3" sets columns; within () set which figures go in which cells
       widths=c(1,1), heights=c(1,2)) # sets ratios of width and height for each cell

plot(mex2, border="grey", axes=F)
plot(mexnbnet2010pan, coordspan, col="blue", pch=".", add=T)
text(coordspan, label=mexpan$ST_ID, cex=.75)
title("PAN 2010") 

plot(mex2, border="grey", axes=F)
plot(mexnbnet2010prd, coordsprd, col="darkgoldenrod1", pch=".", add=T)
text(coordsprd, label=mexprd$ST_ID, cex=.75)
title("PRD 2010")

plot(mex2, border="grey", axes=F)
plot(mexnbnet2010pri, coordspri, col="red2", pch=".", add=T)
text(coordspri, label=mexpri$ST_ID, cex=.75)
title("PRI 2010")

par(mar=c(0,0,1,0))
plot(mex2, border="grey", axes=F)
plot(mexnbnet2010pri, coordspri, col="red2", pch=".", add=T)
text(coordspri, label=mexpri$ST_ID, cex=.75)
plot(mexnbnet2010pan, coordspan, col="blue", pch=".", add=T)
text(coordspan, label=mexpan$ST_ID, cex=.75)
plot(mexnbnet2010prd, coordsprd, col="darkgoldenrod1", pch=".", add=T)
text(coordsprd, label=mexprd$ST_ID, cex=.75)
title("Party Networks (2010)")
par(mar=c(5, 4, 4, 2)+0.1)
dev.off()

# end